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ABSTRACT 

We present results of numerical simulations carried out with a two-dimensional radiation hydrodynamics 
code in order to study the impact of massive stars on their surrounding interstellar medium. This first paper 
deals with the evolution of the circumstellar gas around an isolated 60 M . star. The interaction of the photo- 
ionized H n region with the stellar wind bubble forms a variety of interesting structures like shells, clouds, 
fingers, and spokes. These results demonstrate that complex structures found in H ii regions are not 
necessarily relics from the time before the gas became ionized but may result from dynamical processes dur- 
ing the course of the H ii region evolution. We have also analyzed the transfer and deposit of the stellar wind 
and radiation energy into the circumstellar medium until the star explodes as a supernova. Although the total 
mechanical wind energy supplied by the star is negligible compared to the accumulated energy of the Lyman 
continuum photons, the kinetic energy imparted to the circumstellar gas over the star’s lifetime is 4 times 
higher than for a comparable windless simulation. Furthermore, the thermal energy of warm photoionized 
gas is lower by some 55%. Our results document the necessity to consider both ionizing radiation and stellar 
winds for an appropriate description of the interaction of OB stars with their circumstellar environment. 
Subject headings: galaxies: evolution — H ii regions — hydrodynamics — instabilities — 

ISM: bubbles — ISM: structure 


1. INTRODUCTION 

Massive stars play an important role in the evolutionary 
history of galaxies. They are the primary source of metals 
and dominate the turbulent energy input into the interstellar 
medium (ISM) by stellar winds, radiation, and supernova 
explosions. The radiation field of a massive star first 
dissociates the ambient molecular gas and forms a so-called 
photodissociation region of neutral hydrogen. Subsequently, 
the Lyman continuum photons of the star ionize the H i gas 
and produce an H ii region that expands into the neutral 
ambient medium. A fast stellar wind creates shocks that form 
a so-called stellar wind bubble (SWB) filled with hot plasma, 
which expands into the H ii region. Finally, the star explodes 
as a supernova of Type II (SN II), creating a supernova rem- 
nant (SNR) that sweeps up the ambient medium. The 
SWBs and SNRs of neighboring stars can overlap and form 
a superbubble with a diameter of order 1 kpc. 

This paper pursues two major goals. First, we want to 
examine the combined influence of wind and ionizing radia- 
tion on the dynamical evolution of circumstellar matter 
around massive stars; i.e., we are interested in the inter- 
action processes between the photoionized H ii region and 
the SWB that evolves into the ionized gas. The second goal 
is to improve our knowledge of the energy transfer efficiency 
between massive stars and the ISM: How and to what frac- 
tion is the energy of stellar radiation and the stellar wind 
converted into kinetic, thermal, and ionization energy of the 
ISM? How does the formation of the SWB influence the 


transfer of stellar radiation? To what extent does all this 
depend on the evolutionary state of the star? 

To investigate these effects, we perform numerical two- 
dimensional radiation hydrodynamic simulations of the 
interaction between an isolated massive star and its sur- 
rounding ISM via stellar hydrogen-ionizing photons and a 
stellar wind. We calculate the hydrodynamical evolution of 
the circumstellar gas coupled with radiation transfer, time- 
dependent ionization of hydrogen, and a realistic descrip- 
tion of cooling. The stellar mass-loss rate, the terminal 
velocity of the wind, the effective temperature, and the lumi- 
nosity of the star are specified as time-dependent boundary 
conditions. We examine the evolution of the circumstellar 
material starting from the main-sequence (MS) phase of the 
star until it explodes as a supernova. We do not consider the 
SNR formation. In this paper we present results of these cal- 
culations for a star with an initial mass of 60 M s , extending 
the work of Garcia-Segura, Mac Low, & Langer (1996b, 
hereafter GML1) to a more precise analysis of the energetic 
aspect during the whole evolution of the star. 

The remainder of this paper is struct ured as follows: In § 2 
we briefly review theoretical and observational studies of 
H ii regions and SWBs around massive stars for later com- 
parison. Our numerical procedure, initial conditions, and 
time-dependent boundary conditions are described in § 3. In 
§ 4 we present the numerical results of our calculations and 
discuss them in the context of previous analytical 
approaches, other numerical investigations, and observa- 
tions. Finally, in § 5 we summarize our main results and 
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draw some conclusions from the work presented in this 
paper. 


2. THEORETICAL CONCEPTS AND 
OBSERVATIONAL CONSTRAINTS 


2.1. H ii Regions 

H ii regions are observed in various sizes and shapes. 
Depending on their geometrical extent, they are labeled as 
ultracompact, compact, extended, or giant. The ultracom- 
pact (linear size <0.1 pc) and partly also the compact (0.1- 
0.3 pc) H ii regions are still deeply embedded in their birth- 
places, mature molecular clouds, and as a result of dust 
obscuration generally are only observable in the IR and 
radio wavelengths. On the other hand, extended H ii regions 
with sizes up to several parsecs and giant H ii regions that 
are composed of individual H ii regions and have sizes of a 
few 100 pc are often easily discernible in the optical by their 
bright Balmer and forbidden metal lines. Large holes and 
shells make it obvious that the giant H ii regions are 
powered by not only stellar radiation but also stellar winds 
and supernova explosions in the underlying OB star cluster 
(see, e.g., Yang et al. 1996). 

Within the framework of the simplest theoretical approach 
an O star suddenly “turns on” in a constant-density, 
motionless medium and begins to ionize its surroundings. 
The reaction of the medium to the stellar photons is well 
known and has been described in detail in standard textbooks 
(see, e.g., Spitzer 1978). For comparison purposes with our 
numerical results we are particularly interested in the time 
evolution of the location r, of the ionization front: 


n = 

ro 


3/4 


(Kn'o 

3^Lyc 


T + O 


,7/4 x 4/7 


l>/3 


WoP {2) (T)\ 


( 1 ) 

( 2 ) 


where c s% n is the isothermal sound speed of the ionized gas, 
t\) is the initial Stromgren radius, L Lyc is the photon lumi- 
nosity in the Lyman continuum, n 0 is the hydrogen number 
density in the neutral ambient medium, /3 (2) (T) is the coeffi- 
cient for recombinations of hydrogen into levels 2 and 
higher, and r is the age of the star. According to Lasker 
( 1 967), the kinetic energy of the expanding swept-up shell is 


E k 


2 3/2 r /7 3/4 , 7/4 \ 6/7 3/2i 

t + r 0 ' ) - r 0 ' J , (3) 


where tn H is the mass of a hydrogen atom. Because this 
swept-up shell can be expected to remain cool and neutral as 
a result of strong radiative cooling, we can also estimate the 
ionization energy stored in the H ii region: 


j-* 4 (1 5/2 . 7/2 \ 6/7 

Ei = f 7I7Z0 U c sM r o T + r o ) Xo • 


(4) 


Here we have used equation ( 1 ) for the time- 
dependent radius of the ionization front and \o is the 
ionization potential of hydrogen in the ground state. 

The same Ansatz can then be used to estimate the thermal 
energy of warm gas in the H ii region: 

E, = 4nn 0 (jCv.ii^V + r 7 0 /2 ) (l/1 k B Tn , (5) 

where we will insert T u = 8000 K as an appropriate 
approximation to the temperature in the H ii region: k% is 
the Boltzmann constant. 


2.2. Theory of Stellar Wind Bubbles 

First considerations about the inner structure of SWBs go 
back to Pikefner (1968), Avedisova (1972), Dyson & 
de Vries (1972), and Dyson (1973), but the works of Castor, 
McCray, & Weaver (1975) and Weaver et al. (1977) set a 
milestone in the approach of understanding SWBs. They 
presented a fairly complete picture of the structure and evo- 
lution of SWBs together with a set of equations that 
describes the evolution under the simplifying assumptions 
of a point source of a constant and spherically symmetric 
strong wind that interacts with a homogeneous ambient 
ISM. The global structure that arises from such a wind-ISM 
interaction is depicted in Figure 1 together with the expand- 
ing H ii region into which the SWB evolves. It consists of a 
free-flowing supersonic wind that is heated to about 1 0 6 — 10 8 
K when it passes the inner reverse shock at r s \. The pressure 
of this hot rarefied gas that normally fills most of the volume 
of the SWB is typically much higher than the pressure in the 
photoionized ambient medium. As a consequence, the hot 
gas bubble expands into the H ii region producing a forward 
shock at r S 2 that sweeps up the gas from the H ii region in a 
shell that is separated from the hot bubble interior by a 
contact discontinuity at r c . 

Three phases in the evolution of such a bubble can be dis- 
tinguished. The first is the adiabatic phase that lasts until 
the shock speed v s 2 drops below 200 km s -1 (Falle 1975), 
typically a few times 10 to a few times 10 3 yr, depending on 
the mechanical luminosity of the wind and the ambient den- 
sity. The bubble is expanding so fast that radiative cooling 
does not play a significant role for the dynamical behavior. 
The second stage of evolution is characterized by strong 
cooling in the shell of swept-up material (between r c and 



Fig. 1. — Schematic structure of an SWB: r s \ marks the position of 
the reverse shock, r c the contact discontinuity, r v2 the forward shock of the 
stellar wind bubble, r, the ionization front, and r v3 the forward shock of 
the H ii region expansion. 
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r s2 ), allowing it to be compressed into a geometrically thin 
[r c (t) = r s2 (t), where t denotes time], dense shell, whereas 
for the hot bubble (between r vl and r c ) cooling is still negli- 
gible. This phase lasts much longer than the first, so-called 
adiabatic phase. Thermal conduction from the hot bubble 
interior to the collapsed shell may become important and 
may modify the structure of the bubble in this stage. Weaver 
et al. (1977) assume that an equilibrium between the con- 
ductive energy flux and the mechanical energy flux due to 
evaporation of shell mass in the reverse direction is estab- 
lished. Under the assumptions that in the hot bubble 
between r s . and r, = r s2 the pressure is everywhere the same 
and that the thermal energy contained within this region is 
much higher than the kinetic energy, Weaver et al. (1977) 
derived the following equations for the temporal depend- 
ence of the radius r s2 (t) and the pressure in the hot bubble 
P h (t) in this stage: 


r S 2(t) = r e (t) = l!/ 5 Po 


1 / 5 - 1 / 53/5 


Pb{t) = 


(3850tt) 2/5 ' 


L^pfr^ . 


( 6 ) 

( 7 ) 


Here p 0 = is the density of the ambient medium and 
L w = \ M w vi is the stellar wind luminosity, with M w and v w 
being the mass-loss rate of the star and the terminal velocity 
of the stellar wind, respectively. 

It is also possible to obtain an analytical result for the 
energy that is deposited in the ISM by SWBs. We neglect 
the first (fully adiabatic) stage of bubble evolution because 
it is very short. Cooling in the hot bubble strongly depends 
on the efficiency of heat conduction or how efficiently heat 
conduction is suppressed by magnetic fields. We will assume 
for the analytical theory that there is no heat conduction 
and thus cooling in the hot bubble is unimportant during 
the lifetime of the star. The hot bubble acts like a piston on 
the ambient medium, and analogous to the case of an H n 
region it can be shown that the PdV work is equally distrib- 
uted on kinetic energy of shell motion and thermal energy in 
the shell (which gets immediately lost as a result of cooling). 
If we furthermore assume that the kinetic energy of the 
stellar wind is completely transformed into thermal energy 
at the reverse shock (the strong-jump conditions give 15/16, 
but dumpiness of the stellar wind, which we completely 
neglect, might reduce this value), we get for the kinetic 
energy 


r v^t) 

E k =\ / P h dV a . (8) 

Jo 

V S 2 is the volume inside r s2 . Using equation (6) for r s2 and 
equation (7) for the pressure in the hot bubble yields 

E k =±L w r. (9) 

The thermal energy of the hot gas is thus 

fVsiir) 

E t = L k t — / P b dVa=±L K T. (10) 
Jo 

Besides heat conduction and cooling in the hot bubble, we 
have also neglected the fact that part of the thermal energy 
might be used for collisional ionization of gas. Thus, the 
results for E k and E t are upper limits. 


Although the analytical and semianalytical solutions for 
the evolution of SWBs have been improved in recent years 
(Koo & McKee 1992a, 1992b; Garcia-Segura & Mac Low 
1995; Pittard, Dyson, & Hartquist 2001a; Pittard, 
Hartquist, & Dyson 2001b), a variety of physical effects 
remain to be included in order to achieve a better agreement 
of models and observations. For example, the discrepancy 
between models and observations with regard to the evolu- 
tion of the hot phase in bubbles has recently been reviewed 
by Mac Low (2000; see also Chu 2000). It has become 
evident that the stellar parameters, such as the mass-loss 
rate, the terminal velocity, the effective temperature, and the 
luminosity of the star, that drive the evolution of the circum- 
stellar matter vary strongly over time. Because most pre- 
vious studies dealt with the evolution of either H n regions 
or SWBs, little is known about the interaction of these two 
structures. 


2.3. Observations of Stellar Wind Bubbles 

Hot gas (10 6 — 10 8 K) is expected to be a crucial indicator 
for the existence of SWBs. However, there have been only 
two successful X-ray observations of SWBs so far: NGC 
6888 (Bochkarev 1988; Wrigge, Wendker, & Wisotzki 1994) 
and S308 (Wrigge 1999). Both are actually Wolf-Rayet 
(W-R) bubbles (i.e., the wind-driving star has already 
entered its W-R stage). The fact that no MS bubble has yet 
been observed in X-rays may be due to their physical prop- 
erties; they are expected to be large, diffuse, and conse- 
quently dim (compared with W-R bubbles). There may also 
be a problem with the theoretical models. More sensitive 
X-ray telescopes of the next generation should be able to 
illuminate this problem. 

The shells of MS bubbles are also difficult to observe in 
the optical because they are large and dim (McKee, Van 
Buren, & Lazareff 1984). Nevertheless, shells around MS 
stars, as well as shells around evolved massive stars (which, 
based on their radius, expansion velocities, and shell mass, 
are considered to originate from the MS phase of these 
stars), have been observed in the optical (Lozinskaya 1982; 
Oey & Massey 1994; Marston 1997), the IR (Marston 
1996), and the H i 21 cm line (Cappa et al. 1996; Cappa & 
Benaglia 1998; Benaglia & Cappa 1999). 

The MS shells are typical results of the interaction of 
massive stars with their surrounding ISM over a time 
period of a few megayears. There are also young objects 
that have influenced their interstellar environment so 
strongly that they can already be identified as star- 
gas interactions. One example that should certainly be 
mentioned here is the so-called Becklin-Neugebauer- 
Kleinmann-Low (BN-KL) nebula in the Orion molecular 
cloud 1 (OMC-l) behind the Orion Nebula. It shows a 
spectacular, finger-like outflow of molecular hydrogen 
(Taylor et al. 1984; Allen & Burton 1993; McCaughrean 
& Mac Low 1997; Salas et al. 1999). As a result of 
extremely high visual extinction (A v ~ 20-50 mag), it can 
only be observed in the infrared or at longer wavelengths. 
Most models for the gas outflow proposed so far assume 
that it is driven by mass ejection from a young star, 
namely, IRc2 and/or the Becklin-Neugebauer (BN) 
object, which are both very close to the projected center 
of the mass outflow. IRc2 has a mass-loss rate 
M„. « 10~ 3 to 10 -4 Af s yr _1 (Downes et al. 1981) and a 
luminosity L p hot ~ (2-10) x 10 4 L & (Genzel & Stutzki 
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1989). The mass-loss rate and luminosity of BN are 
M w « 4 x 10“ 7 M 0 yr~* and L pho t ~ (1-2) x 10 4 L© 
(Scoville et al. 1983). 

Another fairly well studied region where formation of 
stars and their interaction with the ambient medium take 
place in an early phase is the Eagle Nebula (Ml 6). 
Imaging with the WFPC2 on board the Hubble Space 
Telescope (Hester et al. 1996) and observations in the 
infrared with ISOCAM (Pilbratt et al. 1998) and 
NIRSPEC (Levenson et al. 2000) reveal a number of 
interesting insights into the structure and possible forma- 
tion history of this active region. The picture of M16 in 
the optical is dominated by the appearance of the so- 
called elephant trunks, impressive pillars of dense molec- 
ular gas protruding into the H n region. Their surface is 
partly illuminated by the ionizing light of nearby (dis- 
tance to the pillars «2 pc) massive stars so that photo- 
evaporative flows from the surface can be identified on 
the well-resolved pictures (Hester et al. 1996) as striations 
extending normally from the surface of the cloud. The 
flow is driven by the pressure gradient that arises when 
the ultraviolet photons from the massive star heat the 
dense gas at the surface of the trunk to a few thousand 
kelvin, resulting in a pressure that is higher than in the 
rest of the (lower density) H n region. 

An interesting question arises in connection with Ml 6: is 
the elephant trunk structure a remnant of the original 
molecular cloud from which the young stars in the region 
were formed, or has this structure actively been molded by 
the massive stars as a shell swept up by stellar winds and 
folded by hydrodynamieal instabilities? More detailed 
observations and better theoretical models are needed to 
answer this question. 

2.4. Previous Numerical Models of Stellar Wind Bubbles: 
An Overview 

Since effects associated with inhomogeneous ambient 
media or time-dependent stellar parameters, as well as 
the formation of gasdynamical instabilities, are difficult 
to handle within an analytical framework, numerical sim- 
ulations have become a powerful tool to study the evolu- 
tion of H ii regions and SWBs. First, two-dimensional 
numerical calculations of SWBs have been performed by 
Rozyczka (1985) and Rozyczka & Tenorio-Tagle (1985a, 
1985b) in a series of papers. For the expansion of the 
bubble into a homogeneous medium the authors 
described the formation of clumps in the thin shell, but 
as a result of the lack of resolution, they were unable to 
completely resolve the dynamics of the shell and to com- 
pare it with theoretical predictions. The authors also 
studied how a change of the wind luminosity (an increase 
between two fixed levels in a certain time) influences the 
shell stability and how SWBs break out from a plane- 
parallel stratified disk. The shell fragmentation due to an 
increase of the wind luminosity has also been studied by 
Stone, Xu, & Mundy (1995) and has been compared with 
the observation of “bullets” in the vicinity of young 
stars. 

GML1 and Garcia- Segura, Langer, & Mac Low (1996a) 
focused on the history of the stellar mass loss and the 
implications for the formation of W-R bubbles. They cal- 
culated the evolution of the SWB in one dimension until 
the star enters the LBV phase (for the model with an initial 


stellar mass of 60 M©) or until the star leaves the red super- 
giant (RSG) stage (for the 35 M© model). The resulting 
one-dimensional profiles of circumstellar density, pressure, 
and velocity were used to set up the ambient conditions for 
two-dimensional calculations of the following stage when 
the W-R wind interacts with the surrounding structure. 
The authors found that in both cases the slowly expanding 
shells originating from prior mass-loss phases are heavily 
eroded when they are overtaken by the faster W-R wind. 
This is also the phase in which the resulting nebula shell is 
considered to be most visible as a result of the high emis- 
sion measure in the high-density filaments resulting from 
the collision. The nebula NGC 6888 with its filamentary 
structure would fit in this picture as being driven by a W-R 
star that has undergone MS and RSG evolution, and the 
shell swept up by the W-R wind is currently colliding with 
the RSG wind shell. 

The approach of Garcia-Segura & Franco (1996) is 
slightly different because they were mainly concerned with 
the very early evolution of H ii regions. They presented two- 
dimensional gasdynamical calculations for the evolution of 
H ii regions in constant and power-law density profiles and 
showed the effect of radiative cooling on the thickness of the 
swept-up shell and therefore also on the development of 
instabilities. They found that the ionization front can rein- 
force the growth of the thin-shell instability, and in a power- 
law density falloff with exponent equal to 2 it can even lead 
to violent shell disruption. 

Brighenti & D’Ercole (1997) investigated the formation of 
W-R shells resulting from fast W-R winds evolving into slow 
anisotropic RSG winds. They found that the shell swept up 
by the W-R wind becomes Rayleigh-Taylor unstable and 
fragments even before it hits the shell that had previously 
been piled up by the RSG wind. X-ray maps computed from 
the numerical models show that the anisotropy of the RSG 
wind leads to a “ two-lobe ” X-ray morphology that is quali- 
tatively in agreement with the observation of NGC 6888 
(Wrigge et al. 1994) but cannot reproduce the small filling 
factor of the X-ray-emitting gas. 

Frank, Ryu, & Davidson (1998) tried a different 
approach: they modeled the interaction of an anisotropic 
fast wind with an isotropic slow wind to explain the mor- 
phology of LBV bubbles. They found that anisotropic fast 
winds can indeed produce strongly bipolar outflows without 
assuming that the fast wind collides with a slowly expanding 
disk or torus as previously postulated. 

Strickland & Stevens (1998) calculated simple MS 
bubbles blown by a wind with constant mass-loss rate and 
terminal velocity into a uniform ISM. Their main goal was 
to calculate synthetic X-ray spectra from these numerical 
models as they would be observed with the ROSA T satellite, 
in order to analyze these spectra according to standard pro- 
cedures and to compare the inferred properties with the 
original numerical models. Surprisingly (or not), they found 
that the inferred properties of the bubble can considerably 
deviate from the “real” properties of the model bubble. 
This implies that detailed X-ray emission models are 
necessary, instead of simple one- or two-temperature spec- 
tral fits, in order to derive reliable properties of bubbles and 
superbubbles. 

In Table 1 we compare the numerical hydrodynamic 
models of SWBs/H ii regions quoted above with respect to 
the included physics, the covered range of parameters, and 
some technical properties. 
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TABLE 1 

Hydrodynamical Models of SWBs/H ii Regions around Single Early-Type Stars 


Linear Size Resolution 3 

References Dimensions Geometry (pc) (xlO -3 ) 


1 2 Cylindrical 7.8-(3.2 x 10 3 ) 5 

2 2 Spherical 0,1 5 

3 1/2 Spherical 70-45/19-3 2.5/2.5-1.25 

4 2 Spherical 0.1-2 5 

5 2 Spherical 6.6/1 1.3 

6 2 Cylindrical 0.2 3.9 

7 2 Cylindrical 6.5 2.5 

8 2 Cylindrical 64 8.0-0.125 


Duration 

References (Myr) Ionization Heating Cooling 


1 0.02-(1.5 x 10 -5 ) Thermal Mechanical CIE b 

2 10 -3 No Mechanical Isothermal EOS c 

3 4-5/(28-7.2) x 10 -3 No Mechanical CIE -F cutoff 

4 0.02-0.33 PIE d Mechanical -I- PIE CIE + cutoff 

5 0.023 No Mechanical CIE 

6 8x 10 -4 No Mechanical CIE 

7 0.015 No Mechanical CIE + cutoff 

8 4 Time dependent 6 Mechanical + radiative Explicit or CIE f 


References Wind Asymmetry L w (t) L Lyc (/) «o( r ) 


1 No Fixed/two-level 0 Constant/composite 

2 No Two-level 0 Constant 

3 No Variable Variable Constant/from one dimension 

4 No 0 Fixed Constant/power law 

5 Yes Two-level with transition 0 Pre-W-R 

6 Yes Three-level 0 Constant 

7 No Constant 0 Constant 

8 No Variable Variable Constant 


a In units of linear size. 

b Cooling function based on the assumption of collisional ionization equilibrium. 
c Equation of state. 

d Photoionization equilibrium calculated assuming that the gas is fully ionized inside the H n region. 
e Thermal and radiative ionization of hydrogen. 

f Explicit calculation of important cooling processes for lower temperatures and CIE for high temperatures. 

References. — (l)Rozyczka 1985; Rozyczka & Tenorio-Tagle 1985a, 1985b. (2) Stone etal. 1995. (3) GML1; Garcia-Segura et al. 1996a. 
(4) Garcia-Segura & Franco 1996. (5) Brighenti & D'Ercole 1997. (6) Frank et al. 1998. (7) Strickland & Stevens 1998. (8) This paper. 


3. NUMERICAL METHOD 
3.1. The Radiation Hydrodynamics Scheme 

The numerical calculations presented in this paper have 
been performed using a two-dimensional radiation hydro- 
dynamics code with a fairly sophisticated treatment of time- 
dependent ionization, radiation transfer, and heating and 
cooling processes. Most features of the code and results of 
test calculations have already been described in Yorke & 
Kaisig (1995) and Yorke & Welz (1996). Therefore, here we 
only give a short summary of the code’s main properties. 

The hydrodynamical equations for the conservation of 
mass, momentum, and energy in an ideal, inviscid, single 
fluid have been formulated in cylindrical coordinates (r, z), 
assuming axial symmetry around the z-axis. This set of 
equations is solved numerically on an Eulerian grid in the 
quadrant with r > 0 and z > 0; i.e., we additionally assume 
mirror symmetry with respect to the equatorial plane. The 
differencing scheme used to discretize the equations is 
second-order accurate in space. Because of operator split- 
ting, the accuracy in time is greater than first order. The 
advection scheme of van Leer (1977) is employed. The star 


that influences the gas is located in the center of the coordi- 
nate system at r = z = 0. We use square grids with constant 
mesh size and insert multiply nested grids in the corner at 
r = z = 0 to enhance the spatial resolution close to the star. 
All the nested grids are self-similar, and the linear spatial 
resolution is improved by a factor of 2 for each level of 
nesting. The basic code is explicit; therefore, the Courant- 
Friedrichs-Lewy (CFL) condition determines the maximum 
time step except when changes in the stellar parameters take 
place on even smaller timescales. Von Neumann-Richtmyer 
artificial viscosity is used for the treatment of shocks. We 
use the equation of state for an ideal gas (in our case a mix- 
ture of molecular, neutral, and ionized hydrogen with elec- 
trons). The formation and dissociation of molecular 
hydrogen is not yet considered in the energy balance 
because it is assumed that continuum radiation below the 
Lyman threshold dissociates the molecules before the gas is 
heated up to temperatures where thermal dissociation 
occurs. 

Our integration order starts with the finest grid. The size 
of the time step satisfies the CFL condition for the finest 
grid, as well as for all coarser grids (normalized by the 


No. 2, 2003 


MASSIVE STARS AND ENERGY BALANCE OF ISM 


893 


appropriate power of 2). After two time steps have been 
done on the finest grid and the solution has been advanced 
to / + 8t\ -F 6/2, one time step 6t[ -f St 2 is done on the next 
coarser grid to advance the solution there to the same time. 
Because of our prior careful choice of and 6/ 2 , their sum 
obeys the CFL condition on the next coarser grid. After the 
solutions on the two grid levels have been advanced to the 
same point in time, all coarse grid values that have under- 
lying fine grid values (not from the boundary) are replaced 
by weighted averages from the fine grid. The outer boun- 
dary conditions of the fine grid are interpolated from the 
corresponding grid values of the coarse grid. After this cycle 
is repeated (taking into account the new CFL time step con- 
straints), the first integration on the next coarser (third- 
level) grid can be performed; third-level values are partly 
replaced by second-level values, and second-level outer 
boundaries are imposed from third-level values, and so on. 
We use this recursive scheme for seven grid levels, which 
means that ^^ =0 2 / = (2 7 - 1) = 127 single integrations 
have to be done before one time step on the coarsest grid is 
completed. 

The ionization structure of hydrogen is calculated for 
each hydrodynamical time step by including the effects of 
photoionization, collisional ionization, and spontaneous 
recombination. The absorption of the stellar point-source 
Lyman continuum by neutral hydrogen is calculated along 
radial lines of sight. We use N r + N- - 1 rays to ensure that 
every grid cell is traversed by at least one ray, where N r and 
N~ are the numbers of grid cells in the r- and z-direction, 
respectively, without ghost cells. The photoionizing radia- 
tion field is assumed to have a blackbody spectrum at the 
effective temperature of the star. We use the gray approxi- 
mation for the transfer of the stellar photons; i.e., we calcu- 
late the absorption cross section and excess energy only for 
the mean energy of photons above the Lyman threshold 
and, therefore, do not account for selective absorption proc- 
esses. To account for the diffuse Lyman continuum that 
originates from electron recombinations directly into the 
ground state of hydrogen, we simply use the on-the-spot 
approximation. 

Radiative heating and cooling are considered as source 
and sink terms in the energy equation and are calculated in 
an additional substep. Photoionization of hydrogen is 
responsible for gas heating. The treatment of cooling 
depends on the gas temperature. Below 1 5,000 K the contri- 
butions by the following processes are explicitly calculated: 
bremsstrahlung, collisional ionization of hydrogen, thermal 
energy lost in hydrogen recombinations, collisional excita- 
tion of neutral hydrogen followed by Lya emission, and col- 
lisional excitation of the low-lying { D terms of N + and 0 + + 
and the 2 D term of 0 + . Simple assumptions have been made 
about the degree of ionization of nitrogen and oxygen: 
because of the similar ionization potentials and comparable 
recombination coefficients of H°, N°, and O 0 , we set the 
ratios of ionized to neutral species for nitrogen and oxygen 
equal to the respective number calculated for hydrogen, i.e., 
Pn + /pn = X and (p Q + 4- po ++ )/po = x, where X is the 
hydrogen ionization fraction. It is furthermore assumed 
that the oxygen ions are equally distributed on the ioniza- 
tion stages 0 + and 0++, i.e., p Q + / po = p Q ++ /po = X /2. 

The chemical composition is solar everywhere in our 
computational domain. Helium is not yet considered 
because it is much less abundant than hydrogen and has no 
low-lying energy levels. For temperatures above 10 5 K we 


use the interstellar cooling function from Sarazin & White 
(1987) with the correction given by Soker (1990). This cool- 
ing function is based on the assumption of collisional ion- 
ization equilibrium. In the intermediate temperature range 
between 15,000 and 10 5 K we calculate both values and 
use a weighted average. We have not included dust in the 
calculations presented here. 

3.2. Initial Conditions 

We start all our calculations with the turn-on of the zero- 
age main-sequence (ZAMS) stellar radiation field and stel- 
lar wind in a homogeneous and quiescent ambient medium. 
This is definitely a gross oversimplification of what can be 
expected for the structure of the circumstellar gas shortly 
after the star has been formed at the end of the pre-MS 
phase. We have two reasons for our approach. The first is 
that, in order to understand what happens in the circumstel- 
lar medium, it is important to start with a well-defined and 
simple initial configuration. As we will show below, even 
with this initially homogeneous ambient density and tem- 
perature distribution, the interaction of the SWB with the 
ionizing radiation field produces a variety of interesting 
morphological structures whose formation processes have 
to be understood in detail before studies with a more realis- 
tic initial setup can be performed. The second reason is that 
we want to keep our simulations comparable with those of 
other authors who use different algorithms or incorporate 
different physical effects in their codes. We choose no = 20 
cm -3 and T 0 = 200 K for our quiescent ambient medium. 
The density is perhaps high for a star that has already left its 
parental molecular cloud, but we want to keep our results 
comparable to those of GML1, who also used no = 20 cm" 3 
in order to reduce computational expenses. The tempera- 
ture T 0 = 200 K in the ambient medium was chosen to yield 
a thermal pressure “typical” for the ISM at the solar 
galactocentric distance. 

We use a spherical one-dimensional solution after 1000 yr 
as the initial model for the two-dimensional simulations in 
order to prevent boundary effects from influencing the for- 
mation of the very small wind bubble while the development 
of the self-similar structure is not completed. We switch 
from one to two dimensions before the swept-up gas 
collapses into a thin shell that is subsequently subject to 
instabilities. 

3.3. Boundary Conditions 

In our models the stellar parameters for mass-loss rate 
(M w ), terminal velocity of the wind (u lv ), effective tempera- 
ture ( T e ff), and photon luminosity in the Lyman continuum 
(L Lyc ) are time-dependent boundary conditions that drive 
and govern the evolution of the circumstellar gas. For the 
60 M e star we adopt the stellar parameters given by GML1 . 
They were obtained from stellar evolutionary models and 
from observations and are shown in Figure 2. The 60 M© 
star is supposed to undergo the following evolution: MS O 
star — » H-rich WN star — > P Cygni-type LBV — ► H-poor 
WN star — > H-free WN star WC star — ► SN. 

According to the mass-loss rate and terminal velocity at 
time /, appropriate values for density and velocity of the gas 
are set within the “ wind generator region ” on the finest grid 
(a small sphere around the center of the coordinate system 
where the star is located). The radius of the “ wind generator 
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Fig. 2. — Time-dependent stellar parameters used as boundary conditions for the calculation of the 60 M s model: mass-loss rate (top left), terminal velocity 
of the wind (top right), effective temperature (bottom left), and luminosity (bottom right). All parameters are adopted from GML1 . 


region” is 3.5 x 10 17 cm for the calculations shown in this 
paper. 

Reflecting boundary conditions are used at the r- and 
z-axes. The outer boundaries of the nested grids were taken 
from the corresponding values on the next coarser grid. The 
outer boundaries of the outermost (coarsest) grid are semi- 
permeable; i.e., mass, energy, and momentum can flow out 
of the computational domain, but not in. Of course, outflow 
would conflict with our intention to carefully take stock of 
the energetic processes in our calculations. Thus, we prevent 
outflow by choosing the coarsest grid large enough that 
neither the ionization front nor moving material can reach 
the outer boundary before the star dies and the calculation 
stops. 

3.4. Geometry and Resolution 

The size of the computational domain is r max = z max = 64 
pc. We use six nested grids within the coarsest grid, resulting 
in seven grid levels. The simulation has been performed with 
125 x 125 cells on each grid level (excluding ghost cells). 
The linear resolution ranges from 8 x 10 -3 pc close to the 
star to approximately 0.5 pc in the outermost parts of the 
coarsest grid, which were affected only during the late stages 
of the bubble evolution. For a resolution study we repeat 
the calculations with 61x61 cells and for the first megayear 
of evolution with 253 x 253 cells. 

4. RESULTS AND DISCUSSION 
4.1.^ Resolution Study 

Before we present our results, we briefly discuss their 
validity. It is currently impossible to spatially resolve all 


details of the entire dynamical evolution when modeling an 
SWB in two or more dimensions over such a long time with 
all the physics we have included. Especially the modes of the 
thin-shell instability (Vishniac 1983; Ryu & Vishniac 1988; 
Vishniac & Ryu 1989) are extremely difficult to resolve 
because a sufficient number of grid cells across the thin shell 
are required to follow the tangential flow of material. Mac 
Low & Norman (1993) performed purely hydrodynamical 
simulations of the thin-shell instability and confirmed the 
linear stability analysis of Vishniac (1983). Since they com- 
pletely focused on the hydrodynamical evolution of the thin 
shell, they made no efforts to treat much more complex 
systems with more sophisticated physics. 

To check the impact of this resolution effect on our results 
for the energetic evolution of the SWB, we performed a res- 
olution study. We calculated the same 60 M 0 model with 
three different resolutions, medium (125 cells in each dimen- 
sion), low (61 cells in each dimension), and high (253 cells in 
each dimension). We compare only the first megayear of the 
evolution (because the high-resolution model is too CPU 
intensive). The results are displayed in Figure 3. The change 
of kinetic energy due to different resolutions is systematic 
(at times later than 0.2 Myr) in the sense that a higher reso- 
lution increases the kinetic energy that remains in the sys- 
tem. For most of the time the difference in kinetic energy 
between medium and high resolution (< 0.05 dex) is smaller 
than between low and medium resolution, which is expected 
as the value converges toward the limit. The variation of the 
thermal energy with resolution is only weak; the maximum 
difference between low and high resolution is <0.05 dex at 
any time during the first megayear. The scatter in ionization 
energy for the three different resolutions is also <0.05 dex 
for 0.5 Myr < t < 1.0 Myr and ;g0.1 dex for t < 0.5 Myr. 
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Fig. 3. — Resolution study for the 60 M 0 model. is the total kinetic 

energy of bulk motion in the system, E t the thermal energy, and £, the 
ionization energy (13.6 eV per ionized hydrogen atom). 

In the latter case, the deviation of ionization energy between 
the medium- and high-resolution runs varies but is con- 
siderably smaller than between low and medium resolution, 
indicating that the value of the ionization energy is already 
close to the actual value. 

These results show that the errors in our energetic 
analysis due to resolution effects are within an acceptable 
range, although some details of the evolution are not yet 
completely resolved in our calculations. 

4.2. The Evolution of the 60 M e Model 

Figures 4-16 show the evolution of the circumstellar gas 
around our 60 M e model star. The plots show gray-scale 
coded physical quantities in a quadrant of the (r, z)-plane 
through the star. The size of the displayed area varies with 
time as the expansion proceeds outward. The star is located 
in the center of the coordinate system. We display the mass 
density together with the velocity field, the gas temperature, 
and the degree of hydrogen ionization, respectively, for cer- 
tain stages of evolution. The degree of ionization is not 
shown for stages after 1 Myr, when the ionized region 
spatially almost coincides with the hot gas. To prevent con- 
fusion, velocity arrows are skipped in the free-flowing wind 
and (except for Fig. 4) in the hot bubble. 

Figure 4 shows the initial model for the two-dimensional 
calculation that has been set up from the one-dimensional 
solution after 1000 yr. The fundamental structure of the 
combined SWB/H n region is well described by the analyti- 
cal solutions. The ionization front at a distance of about ^8 
pc from the star separates the ionized gas at about 8000 K 
from the undisturbed medium. The front has not yet 
reached the initial Stromgren radius at about 12.7 pc and is 
still R type. Thus, there is little dynamical response of the 
heated gas at this point. The SWB evolves into the pre- 
ionized and preheated but otherwise undisturbed medium. 
The free-flowing wind is heated by the reverse shock at 
r«0.16 pc to very high temperatures of about 10 8 K. 
Although the density in the hot gas is fairly low (^4 x 10 -25 
g cm -3 ), the thermal pressure is still much higher than in the 
H ii region and the hot bubble expands supersonically (with 
respect to the isothermal sound speed in the H n region of 
«12 km s -1 ) at roughly 200 km s -1 . The photoionized H ii 


region gas becomes shocked by the expanding hot bubble 
and compressed into a shell with a density of more than 
10~ 22 g cm -3 and a temperature of nearly 10 6 K. At this 
early time, the SWB and the photoionized H n region still 
evolve independently. 

Figure 5 depicts the evolution after 2 x 10 4 yr. The photo- 
ionized region has reached its initial radiative equilibrium 
size, and gas at the outer edge of the H ii region starts to 
accelerate outward as a result of the pressure gradient 
between the warm ionized gas and the cold neutral ambient 
medium. The hot bubble has grown to about 2 pc in radius. 
The shell has been decelerated to approximately 60 km s -1 . 
Because of the high density and high temperature, strong 
cooling has led to compression of the swept-up gas into a 
thin shell. The thin-shell instability initially produces some 
dense knots in the shell that imply variations in the optical 
depth along radial lines of sight indicated by the rippling of 
the photoionization front. The initial evolution of the SWB 
(as long as it is independent of the H ii region evolution) is 
very similar to what is described by Strickland & Stevens 
(1998), except for the deviations due to the different ambient 
density and stellar parameters. 

On the basis of Figure 6 (t — 4 x 10 4 yr) and Figure 7 
(/ = 0. 1 Myr) the structuring impact of the interaction 
processes between the SWB and the H n region can be well 
identified: The photoionized H ii region has started to expand 
into the ambient medium at 6-7 km s -1 . The isothermal 
sound speed of the cold neutral gas is about 1 km s -1 . Thus, 
an outer radiative forward shock occurs that sweeps up the 
ambient gas into a shell. The hot gas bubble has grown to 
about 5 pc in radius, and the shell expansion velocity is still 
about 25 km s -1 ( t = 0.1 Myr). The ionization plots nicely 
illustrate the influence of the density knots in the wind bubble 
shell on the evolution of the H n region: The gas clumps cast 
shadows into the H ii region. The optical depth for hydro- 
gen-ionizing photons along lines of sight through the clumps 
is higher than for lines of sight that pass the shell between 
clumps. As a result of the weakening of the radiation field, 
the gas behind the clumps starts to recombine, forming neu- 
tral spokes within the H ii region. The spokes are less promi- 
nent in temperature than in degree of ionization because the 
cooling efficiency breaks down with the disappearance of the 
free electrons and the shadowed regions remain at several 
thousand kelvin for some time. 

One could imagine that the shadowed regions are artifi- 
cial as a result of the use of the on-the-spot approximation 
for the diffuse radiation field that originates from recombi- 
nations directly into the ground state of hydrogen. How- 
ever, a test calculation that used flux-limited diffusion 
instead of the on-the-spot approximation yielded qualita- 
tively the same results with respect to the shadowed regions. 
Only slight morphological differences were obtained as a 
result of the transfer of the diffuse radiation field. 

Considerations about the influence of shadow patterns 
on ionized regions are not completely new: Capriotti 
(1973) already discussed the impact of shadows by shell 
fragments on the formation and appearance of planetary 
nebulae (PNs). Williams (1999) examined the corrugation 
of R-type ionization fronts by density inhomogeneities 
with subsequent formation of dense clumps. Soker (1998) 
studied the formation of compressed tails in the shadow 
zones behind dense clumps that are present in the vicinity 
of the central star of a PN before the ionization starts. In 
addition. Canto et al. (1998) presented an analytical 
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model for the partial ionization of the shadow regions 
behind neutral clumps by the diffuse radiation produced 
in the nebula along with some gasdynamic simul- 
ations. The model is not restricted to PNs and can be 
generalized to photoionized regions. 


Richling & Yorke (2000) considered the effects of UV 
dust scattering and calculated the evolution of disks exter- 
nally illuminated by an O star. Here, also, shadows are cast 
that appear cometary in shape with tails pointing away from 
the ionizing source. Had we considered UV scattering by 
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Fig. 4. — Circumstellar mass density and velocity field (top), temperature 
(middle), and degree of hydrogen ionization (bottom) for the 60 M . model 
at age 10 3 yr (high-resolution run). The velocity arrows in the free-flowing 
wind zone have been omitted to prevent confusion. The star is located in 
the center of the coordinate system. Please note the dilTerent length scales. 


Fig. 5. --Same as Fig. 4, but at age 2 x 10 4 yr. The velocity arrows in the 
free-flowing wind zone and in the hot bubble have been omitted to prevent 
confusion. 
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Fig. 6. — Same as Fig. 5. but at age 4 x 1 0 4 yr 


dust in this investigation, it is very likely that our clumps’ 
shadows would also have had a cometary shape. 

For the combined H n region/SWB described in our work 
we show that the density inhomogeneities needed to cast the 
shadows can be produced by the action of the stellar wind 
even after the H n region has already been formed in a 
completely uniform medium. In our case the growth of the 


ionized fingers of the H n region is therefore not only due to 
the ionization front instability described by Garcia-Segura 
& Franco (1996). It is strongly triggered by the redistribu- 
tion of mass (and thus opacity) by the action of the stellar 
wind shell. 

The real nature of the spokes, clumps, and instabilities is 
three-dimensional rather than two-dimensional, a fact that 
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Fig. 9. — Same as Fig. 5, but at age 0.5 Myr. Because of the bubble 
expansion, a larger volume is shown. 


Fig. 8. — Same as Fig. 5, but at age 0.22 Myr 


Another consequence for the structure of the H ii region 
and the shell of the SWB can be seen in Figure 8 at 0.22 Myr. 
The disturbed structure of the FI n region, especially the for- 
mation of the neutral spokes and the ionized fingers, has 
induced several flow patterns with velocities up to ^20 km 
s -1 . As a result of the pressure gradient between the neutral 
spokes and the ionized gas, the spokes are compressed by the 


certainly should influence their size distribution, shapes, 
and longevity. Three-dimensional simulations, although 
extremely expensive, would be useful to study these struc- 
tures in more detail. These two-dimensional simulations 
should thus be regarded as qualitative indicators that such 
instabilities exist, a feature missing from one-dimensional 
calculations. 
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Fig. 1 1 . — Same as Fig. 1 0. but for the medium-resolution run 


ionized gas. At the same time the expanding SWB compresses 
the H ii region (and its spokes) as a whole. The clumps in the 
stellar wind shell originally responsible for the shadows grow 
in mass and partially merge during the further course of evo- 
lution. Photoevaporation of the stellar wind shell and its 
embedded clumps does not strongly influence the dynamics 


because the gas in these structures is still warm and its photo- 
ionization does not produce large pressure gradients. How- 
ever, where stellar photons directly impinge upon the cold 
outer shell through the optically thinner ionized fingers, 
photoevaporation occurs and leads to significant dynamical 
effects, as evident in Figures 8-16. 
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Fig. 1 5. — Same as Fig. 1 2. but at age 3.4 1 Myr 


In the top panel of Figure 8 and also in some of the other 
density plots are swirls of material visible in the hot bubble. 
These vortices are sometimes generated at the r- and z-axes, 
possibly as an artifact of the reflecting boundary conditions. 
Since the total mass involved in these flow patterns is very 
small, they do not influence the morphological or the 
energetic evolution significantly. 

Proceeding to 0.5 Myr (Fig. 9; note the different scale), we 
see that the entire H ii region has been swept up by the 
SWB, which has grown to approximately 12 pc in radius. 
The photoionized H ii region (it can best be distinguished 
from the hot bubble using the temperature plot) consists of 
3-7 pc long pillars, the remnants of the ionized fingers. The 
stellar wind shell loses its structure, becomes more or less 
dissolved, and merges into the “remnant” H ii region 
because it sweeps up a highly nonuniform medium with an 
expansion speed that becomes comparable to the sound 
speed of the H ii gas. 

Figures 10 and 1 1 compare the results of our high- and 
medium-resolution runs with respect to the morphologi- 
cal structure after 1 Myr. One can easily see that there 
are significant deviations in the structure of the photoion- 
ized region between the two runs. This is not very sur- 
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prising because, as mentioned above, the formation of 
density fluctuations due to the thin-shell instability is very 
sensitive to the spatial resolution, especially as we have 
started with a completely undisturbed ambient medium 
and have not artificially excited any particular unstable 
mode. We have also seen that exactly these clumps in the 
stellar wind shell strongly influence the evolution of the 
H ii region as a result of their opacity for Lyman contin- 
uum radiation. 

Nevertheless, the general structure is very similar for both 
resolutions: The radius of the hot bubble is roughly 17 pc 
(strong deviations from spherical shape appear in both 
cases). The photoionized gas forms a highly irregular shell 
around the hot bubble with a mean thickness of «5 pc, 
bound on the outside by a swept-up shell of shocked 
ambient material. It appears as if the phenomenon of 
“breakout” of the hot gas into the photoionized gas has 
occurred between the spokes in both cases. However, we see 
no evidence for further breakout into the undisturbed con- 
stant density medium. Had we chosen an ambient environ- 
ment of limited size, i.e., included cloud edges of decreasing 
density within the computational grid, we could expect a 
more dramatic demonstration of breakout. 
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Fig. 1 6. — Same as Fig. 1 2. but at age 4.07 Myr 
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We now jump to t = 3.30 Myr and show in Figure 12 the 
evolutionary state of the circumstellar material shortly 
before the star enters the LBV phase with its strongly 
enhanced mass-loss rate. (Please note that we once again 
changed the plot scale to an appropriate size.) The radius of 
the hot bubble has grown to approximately 40 pc. Because 
of the growth of the bubble (and the slight decrease of H- 
ionizing photon luminosity after 2 Myr), the photoionized 
H ii region has more or less collapsed into a skin with «2 pc 
thickness at the inside of the outer shell. The outer shell itself 
is still corrugated, and so is the photoionized gas. Some 
small clumps have become detached from the shell. 

Figure 13 depicts the structure during the LBV eruption. 
Compared with the earlier plots, one can detect a much 
higher density in the free-flowing wind zone. That is due to 
the higher mass-loss rate of the star and the smaller terminal 
velocity. The rest of the structure is largely unchanged 
because the time step difference is small. 

The LBV phase lasts only for about 10,000 yr and is fol- 
lowed by the final W-R phases of the star, which are once 
again associated with wind terminal velocities of a few times 
10 3 km s _1 . GML1 describes what happens next: the W-R 
wind becomes shocked when it hits the much slower LBV 
wind and the hot, rarefied gas accelerates the dense LBV 
wind (Fig. 14). This is now a highly Rayleigh-Taylor unsta- 
ble situation, and the shocked W-R wind breaks through 
the LBV material, fragments it, and blows it into the MS 
bubble (Fig. 15) where it mixes with the hot, rarefied gas, 
enhancing the density and lowering the temperature there 
for some time. 

The last set of pictures (Fig. 16) shows the structure of the 
circumstellar gas at the end of our simulation. This is the 
pre-SN configuration immediately before the star explodes 
as SN II. The radius of the hot bubble has reached nearly 50 
pc, and the photoionized H n region is limited to the illumi- 
nated inner skin of the outer shell of swept-up ambient gas. 

A “reionization” phase (e.g., Beltrametti, Tenorio- 
Tagle, & Yorke 1982; Tenorio-Tagle et al. 1982) cannot be 
found in our model because the dynamical evolution of the 
circumstellar gas at late stages is dominated by the influence 
of the stellar wind rather than the stellar radiation field. 

4.3. The Energy Balance in the Circumstellar Gas 
4.3.1. Numerical Results 

Figure 17 shows the total energy contributions in the cir- 
cumstellar gas of the 60 M 0 model as a function of time. We 
follow the kinetic energy of bulk motion throughout the 
computational domain, the ionization energy (13.6 eV per 
ionized hydrogen atom), and the thermal energy of cold 
(T < 10 3 K), warm (10 3 K < T < 10 5 K), and hot (T > 10 5 
K) gas. As we have seen before, the ionization and heating 
of the H ii region is the fastest process at the beginning of 
the evolution. Thus, the ionization energy reaches 1.1 x 
10 50 ergs and the thermal energy of warm gas around 
1.5 x 10 49 ergs very soon after the stellar turn-on. 

The evolution of the thermal energy of warm gas follows 
that of the ionization energy over the lifetime of the star 
with a shift of 0. 7-0.9 dex. This is due to the fact that photo- 
ionization by the stellar radiation field is responsible for the 
bulk production of ionized gas at typically 8000 K. The con- 
tribution by shock heating is small because the reverse 
shock does not heat much mass and the thermal energy of 
this gas is accounted under “hot gas” rather than “warm 



Fig. 17. — Temporal evolution of kinetic, thermal, and ionization energy 
in the 60 Af@ model. For details see text. 

gas.” Therefore, we will discuss only the ionization energy, 
bearing in mind that the thermal energy of warm gas 
behaves quite similarly. 

An interesting feature is the drop of ionization energy 
between its first maximum and t « 0.2 Myr. This is due to 
the formation of the dense structures (shell, spokes) in the 
H ii region that reduce the recombination time and increase 
the cooling of the gas (see also Fig. 18). Between t « 0.2 and 
1.9 Myr the denser structures dissolve and the pressure 
decay within the hot bubble together with the H ii region 
expansion result in a decrease of density in the H ii region. 
Thus, the recombination time increases and the ionization 
energy stored in the system increases accordingly. The ion- 
ization energy reaches 2.7 x 10 50 ergs and the thermal 
energy of warm gas 3.7 x 10 49 ergs. 

As the star enters the first WN stage at t « 1.7 Myr, the 
photon luminosity in the Lyman continuum begins to 
decrease slightly while the mechanical wind luminosity 
starts to increase (see Fig. 2). The former directly reduces 
the rate of photons available for photoionization, whereas 
the latter enhances the pressure in the hot bubble, thus com- 
pressing the H ii region, which implies higher recombination 
rates. Thus, ionization energy and thermal energy of warm 



Fig. 18. — Ratio of energies in the 60 M 0 model with and without a stellar 
wind. 
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gas decrease by a factor of 2 and reach local minima at 
t ~ 2.6 Myr. They rise again when the photon luminosity in 
the Lyman continuum increases and the mechanical wind 
luminosity decreases until the star reaches the LBV stage. 

Although spectacular in appearance, the impact of the 
LBV phase on the energy balance is only limited as a result 
of a relatively brief period of time. We note here that the 
table of stellar parameters that we use is sufficiently time re- 
solved even during the LBV phase, contrary to what could 
be suspected from Figure 2. For a detailed plot of the stellar 
parameters during the LBV phase see GML1. The ioniza- 
tion energy in the circumstellar gas drops during the LBV 
stage because the effective temperature of the star decreases 
significantly, but immediately after the LBV eruption, dur- 
ing the next WN stage, the stellar UV flux drives the circum- 
stellar ionization energy to a global maximum at 3.4 x 10 50 
ergs. The ionization energy follows the stellar UV flux. The 
final values of ionization energy and thermal energy of 
warm gas (i.e., at the end of our calculation) are 1.0 x 10 50 
and 2.1 x 10 49 ergs, respectively. 

The kinetic energy of bulk motion in the circumstellar gas 
rises from zero at the beginning (because we started with a 
quiescent medium) as more and more circumstellar gas is 
accelerated by the expansion of the H ii region and the 
SWB. After some 0.3 Myr it reaches the same level as the 
thermal energy of warm gas and rises in the same manner. 
When the mechanical luminosity of the stellar wind 
increases in the first W-R phase, it lifts the kinetic energy 
above 10 50 ergs, close to the ionization energy level. The 
ejection of the LBV nebula raises the kinetic energy for a 
short time almost by a factor of 2, but it drops back when 
the LBV material hits the outer shell and the kinetic energy 
is dissipated into thermal energy. Its final value is 1 .4 x 10 50 
ergs. 

The thermal energy of hot gas accounts for the internal 
energy of the rarefied gas in the SWB that is heated to high 
temperatures by the reverse shock during the MS and W-R 
phases of the star. Since there is no hot gas at the very begin- 
ning, the thermal energy of this gas phase rises from 0. The 
behavior of the curve is very similar to that of the kinetic 
energy, but it increases faster during the first WN stage of 
the star. Between 2.3 and 3.1 Myr it is the dominant form of 
energy in our system, and after some slight changes in and 
around the LBV phase it finally decays with the mechanical 
luminosity of the star to 9.7 x 10 49 ergs. 

The thermal energy of cold gas starts already with the 
internal energy of the quiescent ambient medium. The value 
increases smoothly from 3.4 x 10 49 ergs at the beginning to 
7.9 x 10 49 ergs at the end of the calculation because more 
and more ambient gas becomes swept up by the outer shell. 
The weak shock heats it to several hundred kelvin, whereas 
cooling of this neutral gas is very low. 

At the beginning the radiative energy input has the stron- 
gest influence on the energy distribution. The dynamical 
response of the circumstellar gas (the expansion of the H ii 
region and the acceleration of the material by the SWB) 
takes much longer. Thus, the ionization energy dominates 
the circumstellar energy distribution over the first ^2 Myr. 
After the kinetic energy and thermal energy of hot gas have 
caught up with the thermal energy of warm gas, these three 
forms of energy are comparable between 0.3 and 1.8 Myr. 
Later, because of the accumulation of kinetic energy in the 
shell and thermal energy of hot gas in the bubble, these two 
become comparable with the ionization energy. Between 2.2 


Myr and the end of the calculation these three forms of 
energy differ by not more than a factor of 3, although some 
changes occur. The LBV phase of the star induces very rapid 
changes in the relative energy ratios of at most a factor of 3, 
but they are transitory. 

It is interesting to compare the energy stored in the system 
for the two cases: (1) with wind and (2) the pure H n region 
evolution without wind (Fig. 18). Obviously, the total 
kinetic energy in the calculation with wind is always higher 
than in the corresponding calculation without because of 
the added kinetic energy of the SWB shell. At the end of the 
calculation, the kinetic energy in the model with stellar wind 
is 4 times higher than in the model without. 

By contrast, the ionization energy reaches the same level 
of wl.l x 10 50 ergs very soon after the start of the calcula- 
tion in both cases. As we have already seen in Figure 17, the 
ionization energy in the model with wind drops by a factor 
of 2 during the first 0.2 Myr, whereas in the model without 
wind it continues to grow. This supports our explanation 
that inhomogeneities in the H n region (shell, spokes) result 
in shorter recombination times and thus lead to fewer 
photoionized hydrogen atoms in the system. Although there 
is additional energy input by the stellar wind, the ionization 
energy is lower than in the calculation without wind by 
0. 1-0.4 dex for most of the evolution. This, however, does 
not imply that the Ha luminosity is different. The Ha 
luminosity directly mirrors the stellar ionization input. 

Although the thermal energy of warm gas (55% below the 
respective value of the windless model at the end of the cal- 
culation) is similarly affected as the ionization energy, the 
total thermal energy is always higher in the model with wind 
because we have additional energy deposited into the hot 
and into the cold gas phases. 

Whereas the mechanical power of the stellar wind and the 
radiative power of the star are the energy sources in our 
system, the only energy sink term is radiative cooling of the 
gas. All other energetic processes lead to a redistribution of 
forms of energy within the system but not to a change of the 
total energy. Figure 19 shows the total energy loss of the sys- 
tem due to cooling of the gas and distinguishes three differ- 
ent contributions: 66 H ionization loss ” refers to the 13.6 eV 
binding energy per hydrogen recombination into all states 
above the ground state, while “ H thermal loss ” accounts 
for the thermal energy of the electrons lost by the same proc- 
ess; “ other processes ” include all other cooling processes 



Fig. 19. — Energy loss due to cooling in the 60 M e model 
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Fig. 20.— Energy transfer efficiency with respect to the total energy input 
in the 60 M ( ., model without wind. 


Fig. 21 . — Energy transfer efficiency with respect to the total energy input 
in the 60 M 0 model. 


considered, such as collisionally excited line emission and 
bremsstrahlung. One can see from Figure 19 that the shape 
of all the curves follows the Lyman continuum photon lumi- 
nosity of the star (see Fig. 2). The total energy-loss rate dur- 
ing the MS phase is close to 10 39 ergs s _1 , approximately the 
same as the stellar Lyman continuum luminosity. Although 
the ionization energy in the system rises to 10 50 ergs during 
the first few times 10 4 yr, it is only a small fraction of the 
total input energy. 

“FI ionization loss” is the dominant energy sink term 
during the evolution. It accounts for roughly 80% of the 
total emitted power during the MS and first WN stage of the 
star. After the LBV phase, the relative importance of “ other 
processes ” has risen because the LBV ejecta enhanced the 
density in the hot bubble where collisionally excited lines 
and bremsstrahlung dominate the energy emission. 

Figure 20 shows the transfer efficiency into kinetic, ion- 
ization, and thermal energy and their sum for the model 
without wind. We define the transfer efficiency as the cumu- 
lative fraction of the input energy that has been converted 
into a particular form up to the time t = r. In this case we 
have only radiative input energy and all values are related to 
that. 

The transfer efficiency into ionization energy has the 
highest value during the whole calculation. The transfer effi- 
ciencies into ionization and thermal energy strongly drop at 
the beginning of the calculation. This is a consequence of 
the fact that the photoionization equilibrium has been 
quickly established and the ionization energy stays constant 
within an order of magnitude, whereas the transfer effi- 
ciency is related to the accumulated input energy that grows 
continuously in time. The transfer efficiency into ionization 
energy reaches x 10 -3 immediately before the LBV 
phase and drops as a result of the decreasing photon lumi- 
nosity in the final W-R phases by almost a factor of 5. The 
transfer efficiency into thermal energy reaches the level of 
10 -3 before the LBV phase but decreases only by less than a 
factor of 2 during the remaining evolution. The discrepancy 
between the declines of both values can again be explained 
by the fact that cooling in formerly photoionized gas is 
greatly reduced as a result of missing electrons after recom- 
bination occurs. Thus, the recombined neutral gas has lost 
all its ionization energy but can retain a considerable por- 


tion of its thermal energy over a relatively long period. The 
transfer efficiency into kinetic energy has reached a nearly 
constant level of ^3 x 10~ 4 after 1 Myr and remains at 
that level until the end of the calculation; i.e., after the first 
megayear the kinetic energy in the system grows almost 
proportionally to the total input energy. 

Now we check the impact of the stellar wind on the energy 
transfer in the system. Figure 21 shows the same transfer 
efficiencies as Figure 20 but for the model with wind and 
radiation. Whereas for the model without wind the transfer 
efficiency into ionization energy dominated the total energy 
transfer efficiency, for the model with wind this is only true 
during the first ^2.2 Myr. In the first WN phase the transfer 
efficiency into thermal energy dominates for about 1 Myr 
because the photon luminosity in the Lyman continuum 
decreases and the mechanical wind luminosity increases, 
thus enhancing the production of hot gas. Curiously, in the 
LBV phase all three energy transfer efficiencies have the 
same value of about 2 x 10 -3 . As in the case without wind, 
the transfer efficiency into ionization energy drops more 
strongly after the LBV phase than that into thermal energy. 
The drop of the latter one in this model is mostly due to 
energy loss of the hot gas that dominates the thermal energy 
at late stages of the evolution. 

The total energy transfer efficiency reaches approximately 
4 x 10 -3 at the end of the calculation, about twice as much 
as in the model without wind. This is interesting because we 
know from Figure 2 that the mechanical wind luminosity 
integrated over the whole lifetime of the star is almost negli- 
gible compared to the total radiative energy input in the 
Lyman continuum. In other words, although the stellar 
wind itself does not actually inject a considerable amount of 
energy into the circumstellar gas (compared to the stellar 
radiation field), its presence almost doubles the total energy 
that is finally contained in the gas. This difference is due to 
the fact that the thermal energy deposited increased by a 
factor of 2.4 and that the kinetic energy rose by a factor of 4 
compared to the model without wind. The ionization energy 
actually decreased slightly by some 20%-30%. 

To facilitate the discussion in the next section, we summa- 
rize the values of the individual energy components at the 
end of the simulation in Table 2. Besides the kinetic energy 
of bulk motion ( E k ) and the ionization energy of hydrogen 
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TABLE 2 

Energy Components at the End of the Simulation 


E k E( Ej' C old -£/,w arm ^/,hot 

Model (xl0 49 ergs) (xl0 49 ergs) (xl0 49 ergs) (x!0 49 ergs) (xl0 49 ergs) 


Windless 3.5 13 5.5 4.7 0 

With wind 14 10 7.9 2.1 9.7 


Note. — The thermal energy of the cold component, £, <C oid, contains the internal energy of the initially 
unperturbed ambient medium (3.4 x 10 49 ergs) that has to be subtracted whenever the input of thermal 
energy into the system is considered. 


(Ei), we list the thermal energy of the cold, warm, and hot 
gas components (E'r.coid? E t , warm, and E^hot) f° r the two 
models, one with and one without stellar wind. The values 
of the energy transfer efficiency into kinetic energy (e*), ion- 
ization energy (e t ), and thermal energy (£,) at the end of the 
simulation are given in Table 3. The total energy of the 
Lyman continuum radiation emitted by the star is 
£ L yc = 1.07 x 10 53 ergs, and the mechanical energy injected 
into the system by the stellar wind, if considered, amounts 
to E w = 3.31 x 10 51 ergs. 


4.3.2. Comparison with Analytical Results 

We now calculate the values for the kinetic, ionization, 
and thermal energy in the system according to the analytical 
solutions given in § 2. The analytical approach cannot 
handle time-dependent stellar parameters; therefore, we 
simply choose mean values for effective temperature and 
luminosity in the Lyman continuum over the lifetime of our 
model star. With (r e ff) = 5.03 x 10 4 K, {Llvc)=8.33x 
10 38 ergs s -1 , and using ap = 3.37 x 10~ 13 cm 3 s~* as the 
hydrogen recombination coefficient and c Si u = 1.15 x 10 6 
cm s _1 for the isothermal sound speed in the H n region 
(corresponding to T\\ = 8000 K), we obtain the following 
for no = 20 cm -3 after r = 4.065 Myr from equations (3), 
(4), and (5) for the model without wind: 

E k = 4.3 x 10 49 ergs , 

Ej = 5.0 x 10 50 ergs , 

E t = 7.7 x 10 49 ergs . 


The energy transfer efficiency in the analytical approach 
is then defined as 


5 = 


E 

Tilt yC ) ’ 


( 11 ) 


where E can be any of E k , E h and E h depending on the 
efficiency that is calculated. Thus, we obtain 

e k = 4.0 x 10' 4 , 

e t = 4.7 x 10~ 3 , 

e, = 7.2 x 10' 4 . 


TABLE 3 

Energy Transfer Efficiencies at the End 
of the Simulation 


£ k £j £, 

Model (xlO- 4 ) (xlO- 4 ) (xl0“ 4 ) 


Windless 3.3 12 6.4 

With wind 13 9.1 15 


Bearing in mind all the assumptions and approximations 
that have been made to obtain the analytical expression for 
the kinetic energy of bulk motion (see Lasker 1967), the 
deviation of less than 30% in e k is not too bad. This differ- 
ence is not due to the constant effective temperature and 
luminosity that we have chosen to obtain the analytical 
result. Recalculating e k for the first megayear, where the stel- 
lar parameters are about constant, makes the discrepancy 
between the analytical and numerical value of e k even larger. 
The difference is more likely due to temperature deviations 
from 8000 K in the H ii region of the model calculation. 

Comparing the analytical and numerical results for the 
transfer efficiency into ionization energy shows fairly good 
correspondence immediately before the LBV phase. At the 
end of the calculation the numerically determined transfer 
efficiency is about 0.6 dex below the analytical value because 
the Lyman continuum luminosity of the star drops signifi- 
cantly during the final W-R phases of the star and the ion- 
ization energy follows immediately. Of course, the analyti- 
cal solution cannot reproduce this feature. The correlation 
between the drop of Lyman continuum luminosity and the 
thermal energy is weaker because cooling ceases when the 
plasma becomes neutral, as already mentioned above. Thus, 
the correspondence between analytical and numerical trans- 
fer efficiency into thermal energy is much better and the 
deviation at the end of the calculation is less than 15%. 

It is more difficult to compare analytical with numerical 
results in the case of the combined H n region/SWB model 
because we have only the analytical energy transfer effi- 
ciency solutions for the H ii region and SWB separately. We 
simply add up the energy contributions from the H ii region 
and the SWB bearing in mind that this is only a rough 
approximation, which actually neglects the interactions 
between both structures. Because we have not considered 
cooling in the hot bubble, the analytical energy transfer 
rates into kinetic and thermal energy are upper limits. 

We insert the mean value of the mechanical wind lumi- 
nosity, (L w ) = 2.58 x 10 37 ergs s -1 , into equations (9) and 
(10), and together with the results for the pure H ii region 
we obtain 


E k = 9.5 x 10 50 ergs , 

Ei = 5.0 x 10 50 ergs , 

E t = 1.6 x 10 51 ergs . 

To find the energy transfer efficiency, we relate these 
values to the sum of Lyman continuum radiation energy 
and the mechanical wind energy (which, of course, is almost 
negligible), 

£ 

6 = r«Z Lyc ) + (Lj) ’ (12) 
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where E can again be any of E k , E h and E t . We get 
e k — 8.6 x 10 -3 , 

£i = 4.6 x 10" 3 , 

£ t = 1.4 x 10" 2 . 

We have already discussed the fact that the ionization 
energy in the numerical model with wind is lower than in the 
model without wind. This slightly enhances the discrepancy 
between the analytical and the numerical results for the ion- 
ization energy in the case with stellar wind. On the other 
hand, there is no indication that the neglect of the ionization 
energy in the hot bubble is a bad approximation. 

The increase of the thermal energy deposit from the wind- 
less H ii region model to the combined SWB/H n region 
model (factor 2.4) is almost an order of magnitude below the 
analytical upper limit. Since the conversion of thermal energy 
by PdV work is the major source of kinetic energy, the 
increase of the kinetic energy deposit between the models 
without and with wind (factor of 4 at the end of the calcula- 
tion) is also 82% below the analytical upper limit. These find- 
ings are supported by the observation that the radius of the 
bubble at the end of the calculation (^50 pc) is considerably 
smaller than the analytical result of ^68 pc according to 
equation (6). One reason for the lack of thermal energy in the 
numerical model is the variation of the mechanical luminos- 
ity of the star. The LBV wind enhances the density in the 
bubble, which leads to stronger cooling, and the mechanical 
luminosity of the stellar wind generally decreases during the 
last 0.6 Myr. But also during the other stages of evolution, 
cooling in the hot bubble is not completely negligible for the 
energetics of the system. The thermal energy of the hot bub- 
ble calculated analytically at t = 1 Myr is already a factor of 
4.5 larger than the corresponding value in the numerical 
model. Although we have not implemented heat conduction, 
effects like turbulent mixing between hot and cold gas and 
numerical diffusion enhance the cooling rate. The resolution- 
dependent numerical diffusion is not a key factor as we have 
shown in our resolution study. However, the cooling due to 
gas mixing is enhanced because the SWB expands at the 
beginning into the highly nonuniform H n region, a fact that 
cannot be properly handled by the analytical description. In 
general, the analytical results for the upper limits of the trans- 
fer efficiency into thermal and thus kinetic energy of bulk 
motion are much higher than the values from the numerical 
simulations even without heat conduction. 

4.4. Direct Observa t ional Im p l ica lions 

Finally, we briefly discuss some direct observational 
implications of our models. Since for a correct construction 
of intensity maps one would need three-dimensional 
models, the results of our two-dimensional calculations can 
only give a rough estimate of the observable intensities. 
Thus, we have calculated angle-averaged intensity profiles 
for Ha and soft X-ray emission based on the cylindrical 
symmetry in our models. 

The Ha emissivity is calculated according to the table 
(case B) in Osterbrock (1989). We use a r n 0 - 92 fit for the 
temperature dependence and set the emissivity to zero where 
the degree of hydrogen ionization is below 10 -4 . We do not 
account for absorption. 

In Figure 22 we plot the angle-averaged Ha intensity pro- 
files at two different times (t = 0.22 Myr, during which 



Fig. 22. — Angle-averaged Ha intensity for the 60 M s model with and 
without a stellar wind compared at two evolutionary times. 


intense structure formation occurs in the H n region, and 
t = 3.365 Myr, at a late stage after the LBV phase of the 
star). For comparison, we have also plotted the Ha intensity 
profiles of the corresponding models without stellar wind 
for the same times. 

At t — 0.22 Myr the pure H ii region model without 
stellar wind shows the intensity profile of a spherical, emit- 
ting volume that has just started to expand. On the other 
hand, the intensity profile of the model with stellar wind 
shows a dominant peak between r = 8 and 9 pc. This peak 
represents the global intensity maximum at that time and is 
produced by the emission of the dense shell swept up by the 
hot bubble of shocked stellar wind gas. A secondary, much 
smaller peak between r = 13 and 14 pc originates from the 
trapping of the ionization front in the dense outer shell 
fragments piled up by the expanding H n region. 

Comparing these data with the intensity profiles at 
t = 3.365 Myr shows that the further expansion of the 
bubble and H ii region generally lowers the Ha surface 
brightness, as expected for comparable Lyman continuum 
fluxes. In particular, when the SWB completely overtakes 
the H ii region and enlarges the whole structure, the Ha sur- 
face brightness is lower than in the pure H ii region model at 
the same time. 

As a result of the facts that the photoionized region is 
more or less limited to the illuminated inner part of the shell 
and that the Ha emissivity in the hot cavity is very low, the 
intensity profile of the combined SWB/H ii region appears 
slightly limb brightened. The Ha emission from the W-R 
bubble is barely visible above the Ha background from the 
MS bubble. It might be that in our simulation the shell of 
the W-R bubble is too hot (e.g., as a result of the neglect of 
heat conduction) and thus emits more strongly in X-rays 
(see below) than in the optical. 

For the X-ray intensity profiles in the energy range from 
0.5 to 3.0 keV we use the emissivity calculated with the 
Raymond & Smith (1977) program for cosmic chemical 
composition (Allen 1973). The emissivity is set to zero for 
temperatures below 10 5 K; absorption is not considered. In 
Figure 23 we display the angle-averaged soft X-ray intensity 
profiles at t = 1.0, 3.30 (before the LBV phase), and 3.365 
Myr (after the LBV phase). The two snapshots before the 
LBV stage show intensities between 10~ 10 and 10~ 8 ergs s _1 
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Fig. 23. — Angle-averaged soft (0. 5-3.0 keV) X-ray intensity for the 
60 Mq model with stellar wind at selected evolutionary times. 

cm“ 2 sr -1 for lines of sight through the hot bubble. After the 
LBV stage, the soft X-ray intensity close to the center of the 
bubble is strongly enhanced by more than 3 orders of mag- 
nitude to values above 10~ 6 ergs s' 1 cm -2 sr _l when the 
shocked W-R wind hits the LBV ejecta. The soft X-ray 
emission comes mainly from the inside of the shell of the 
swept-up LBV wind, before it breaks apart. 

The W-R nebula RCW 58 around HD 96548 (=WR 40) 
is a possible candidate for an observed counterpart of the 
W-R nebula that forms in our model calculation at 
t — 3.365 Myr. The W-R star is currently of type WN8 and 
there are indications that the star passed the LBV stage 
(GML1; Humphreys 1991). The mean radius of the optical 
nebula is 3.5 (Chu 1982) or 2.5 pc (Arthur, Henney, & 
Dyson 1996) depending on the estimation of the distance to 
HD 96548. Thus, the evolutionary state of the nebula 
should be approximately comparable to the nebula in our 
model. 

Although the strong increase of the surface brightness in 
our model after the LBV phase is basically in agreement 
with the fact that up to now only W-R bubbles (and no MS 
bubbles) have been observed in X-rays, the comparison with 
RCW 58 shows that our model also suffers from the same 
problem that all analytical and numerical models thus far 
have: the X-ray luminosity is much higher than observed. In 
our model, the X-ray luminosity of the W-R bubble in the 
energy range from 0.5 to 3.0 keV is 7 x 10 33 ergs s -1 , 
whereas Moffat et al. (1982) report that no X-ray emission 
from RCW 58 was detected with the Einstein satellite, mean- 
ing that the X-ray luminosity of RCW 58 in this energy 
range must be below 10 32 ergs s -1 (Arthur et al. 1996). 

4.5. The Impact of the Ambient Medium the Stellar 
Parameters , and Unconsidered Physics 

Because of the immense computational efforts, we can 
only show results for one set of ambient medium parameters 
(no = 20 cm~ 3 and To — 200 K). We shall utilize the scaling 
laws for H n regions in ionization/recombination and 
heating/cooling equilibrium with negligible gravity (see 
Table 2 of Yorke et al. 1982) and briefly discuss the changes 
of the results that we expect if other values for the ambient 
medium are used. Separating the evolution of the H n region 
and SWB, we note that if the density in the ambient medium 


is modified by a factor 6 while assuming the same equili- 
brium temperature and same ionizing luminosity, all linear 
scales A_ and timescales f transform according to 
A oc f oc S ~ 2 / 3 , velocities remain unchanged, and mass out- 
flow rates and mechanical luminosity (as well as other rates 
of energy changes) transform as <5 -1 / 3 . Values for the vari- 
ous energies given at a particular time scale as Thus, an 
increase in ambient density results in smaller H ii regions, 
lower energy injection rates, and proportionally lower 
values for the ionization, thermal, and kinetic energy. These 
same scaling laws can be derived in a less rigorous manner 
from our equations (l)-(5). 

Modifying the temperature of the ambient medium does 
not strongly affect the evolution as long as the gas is not ion- 
ized and the pressure in the H n region is much higher than 
the pressure of the ambient medium. 

With respect to the SWB, equation (6) shows that the 
radius of the SWB is proportional to p^ 5 ; i.e., at com- 
parable times it is smaller (larger) by <5*/ 5 in the case of an 
ambient density higher (lower) by a factor 6. In contrast to 
the H ii region, the growth of the thermal and kinetic energy 
of the SWB does not depend on the ambient density, as long 
as the assumptions under which equations (9) and (10) were 
derived remain fulfilled (e.g., thermal pressure in the hot 
bubble is much higher than the thermal pressure in ambient 
medium). 

It is difficult to estimate the impact of a variation of the 
ambient conditions on the highly nonlinear H ii region/ 
SWB interaction effects. Basically, an increase of the density 
in the ambient medium leads to both a smaller SWB and a 
thicker shell of swept-up material. Thus, the ratio of shell 
thickness to SWB diameter grows, which increases the 
unstable length scales and tends to inhibit thin-shell instabil- 
ities with all their associated effects. On the other hand, 
dumpiness in the circumstellar medium may exist anyway, 
abolishing the need to trigger the interaction effects by shell 
instabilities. It is clear that numerical simulations under a 
variety of different ambient conditions are highly desirable 
in order to study these effects in more detail. They will 
become available with increasing computer power. 

It should also be noted here that the stellar parameters 
that we use as time-dependent boundary conditions in our 
simulations are still somewhat uncertain. For example, 
Martins, Schaerer, & Hillier (2002) showed that the determi- 
nation of the effective temperature of O dwarfs based on 
non-LTE line-blanketed atmosphere models including 
stellar winds results in differences to previous calibrations 
up to a few times 10 3 K. The mass-loss rate and its temporal 
variation in the short and violent evolutionary stages like 
the LBV phase are probably even less well known. We have 
chosen the set of stellar parameters from GML1 because it 
is currently state of the art, sufficiently time resolved, and 
we wanted to directly compare our two-dimensional 
simulations to their one-dimensional calculations. 

Nevertheless, deviations of the stellar parameters from 
this set will certainly influence the results of the model 
calculations. Returning to the scaling laws discussed by 
Yorke et al. (1982) and scaling the hydrogen-ionizing 
photon luminosity by / but keeping density and tempera- 
ture constant, we find that lengths and times scale 
according to A oc f ex Z 1 / 3 , velocities remain unchanged, 
rates of energy changes scale like / 2 / 3 , and values for 
energy injected into the ISM at a given time scale like /. 
In a similar manner, a higher stellar wind luminosity will 
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increase the size of the SWB and the energy transferred 
to the shell at a given time. Note that although the ana- 
lytical description of the SWB evolution depends on the 
stellar wind luminosity, i.e., only on the product of mass- 
loss rate and wind terminal velocity squared, this does 
not account for cooling in the hot bubble, which will be 
enhanced when the mass-loss rate rises without a propor- 
tional increase in wind velocity. 

The impact of changes in the temporal evolution of the 
stellar parameters (e.g., a different number of outbursts of 
varying length) on the results of our calculations is difficult 
to estimate. Highly nonlinear effects appear when, e.g., a 
slow shell is overtaken by a fast wind and instabilities arise. 
Interaction processes between the SWB and the H n region 
may also be amplified or damped by variations of the stellar 
parameters. Numerical simulations for different stellar 
evolutionary histories are necessary to study these effects in 
detail. 

Some physical processes are not yet considered in our 
models, and we want to briefly discuss their possible influ- 
ence here. One of these processes is heat conduction, which 
may play a role for the energy transfer in the system. 
Although classical theory and our numerical models with- 
out heat conduction predict that the cooling time for the 
interior of the MS bubble is longer than the lifetime of the 
star, there is indirect evidence that, at least in NGC 6888, 
the MS bubble is cold at the time when the W-R wind blows 
out into the MS bubble (Mac Low 2000). A possible explan- 
ation is mass loading of the MS bubble: heat conduction 
evaporates cold material from the shell or from immersed 
clouds, which enhances the density in the bubble. The 
higher density reduces the cooling time of the bubble. Such 
evaporation from a cloud embedded in hot interstellar gas 
has recently been observed (Chu 2000). We will discuss the 
case of NGC 6888 in a forthcoming paper that describes 
numerical simulations of the circumstellar medium around 
a 35 M 0 star, a model that is more appropriate for 
NGC 6888. 

On the other hand, magnetic fields provide an efficient 
mechanism to reduce the efficiency of heat conduction. In a 
plasma the electrons are constrained to move along mag- 
netic field lines. If the magnetic field lines are aligned with 
the shell or if they are tangled, the effective path lengths for 
fast electrons moving from the hot to the cold gas are much 
longer and heat conduction by electrons is less efficient. Sat- 
uration effects of heat conduction due to electric fields aris- 
ing from charge separation further complicate the situation. 
All these effects are worthy to be studied in detail, but the 
simultaneous inclusion into multidimensional hydrodynam- 
ical simulations remains a computational challenge for the 
future. 

We only allowed for isotropic winds in our models. How- 
ever, line-driving theory for rotating stars suggests that 
mass outflow from the poles should be both more vigorous 
and faster. As already mentioned in § 2.4, Brighenti & 
D’Ercole (1997) and Frank et al. (1998) showed that aniso- 
tropic winds can indeed influence the formation and appear- 
ance of circumstellar shells, producing lobes or bipolar 
outflows. Similar effects may be expected for our models 
when anisotropic winds are employed, but details remain to 
be explored. 

Other physical effects not yet considered in our models 
include variations of the metallicity in the circumstellar gas 
and the photodissociation of molecular hydrogen. In its 


later stages, the star is expected to eject gas that is enriched 
with metals. This enhancement of the metallicity in the 
circumstellar medium leads to stronger cooling and thus 
influences the circumstellar evolution. The dissociation of 
molecular hydrogen by stellar photons with energies below 
the Lyman threshold adds additional (exterior) layers to the 
basic structure shown in Figure 1. The H n region evolves 
into the photodissociated region, which itself expands into 
the molecular gas. Interaction processes between the photo- 
dissociated region and the H n region similar to those 
between the H ii region and the SWB may arise. 

Finally, large-scale density gradients such as those 
encountered at the edge of the clouds will give rise to a 
variety of additional dynamical effects, e.g., champagne 
flows of H ii regions and breakout of hot gas, as well as the 
interactions between these two phenomena. 

5. SUMMARY AND CONCLUSIONS 

Our numerical models for the evolution of circumstellar 
gas around a 60 M s star show that the interaction of the 
stellar wind bubble with the stellar radiation field can 
strongly influence the morphology of the circumstellar 
medium. The rearrangement of circumstellar gas by the 
stellar wind influences the way it reacts to the ionizing radia- 
tion. Density clumps formed in the shell of gas swept up by 
the stellar wind bubble cast shadows into the H ii region. The 
resulting pressure gradients force material into the shadowed 
regions, enhancing their density and forming neutral 
“ spokes,” which subsequently raise the mass of the clumps in 
the shell. The H ii region is extended in directions free of 
clumps. Thus, the formation of these elongated H ii region 
“fingers” in our model not only occurs as a result of the 
ionization front instability described by Garcia-Segura & 
Franco (1996) but is triggered and amplified by the 
redistribution of mass by the action of the stellar wind shell. 

These results also shed light on the open question of 
whether the complex structures that can be found in H ii 
regions are primordial, i.e., relics from the time before the 
gas became ionized, or formed by dynamical processes in 
the course of the H n region evolution. While there are 
strong observational indications that the former plays an 
important role, our results give support to the idea that the 
latter cannot be completely excluded: intense structure 
formation in H ii regions with strong stellar winds can 
occur even if the neutral ambient medium was initially 
homogeneous. 

Nevertheless, these structures are only a temporary 
phenomenon because the extended H n region is eventually 
swept up by the stellar wind shell. If we compare our results 
at t = 3.30 Myr with the one-dimensional solution of 
GML1, we see that, from a morphological point of view, 
there is basically little difference in the overall structure 
except for the appearance of an H ii region at the inner part 
of the remnant shell and a few filaments in the hot bubble. 
Thus, the approach of GML1 to use the one-dimensional 
solution for the estimation of the initial conditions for closer 
studies of the LBV and subsequent W-R stage appears to be 
valid. 

The structure formation induced by the interaction of the 
stellar radiation field with the SWB also has implications for 
the energy balance of the circumstellar gas, mostly via the 
effect that denser gas has a shorter recombination time and 
a higher cooling efficiency. This can be seen in the decrease 
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of ionization and thermal energy of warm photoionized gas 
in the circum stellar medium during the time when the struc- 
ture formation occurs in the first few times 10 5 yr. Further- 
more, the ionization energy lags behind the corresponding 
value in the model without wind by 0. 1-0.4 dex and the ther- 
mal energy of warm gas by 0.2-0. 6 dex for almost the entire 
lifetime of the star. This also implies that the H n regions 
around stars with winds can have a higher emission measure 
than undisturbed ones. 

The ionization energy and the thermal energy of warm 
gas in the circumstellar medium are lower compared to the 
models without wind. Although the total mechanical wind 
energy is negligible compared to the accumulated energy of 
the Lyman continuum photons, the simulation with a stellar 
wind results in a kinetic energy in the circumstellar medium 
that is 4 times higher than in the model without wind. The 
total thermal energy is almost twice as high (or, subtracting 
the initial thermal energy of the background gas, the 
thermal energy that is added to the system is enhanced by a 
factor of 2.4). The energy transfer efficiency of the stellar 
Lyman continuum radiation over long timescales is so low 
because this radiative energy is mostly used to maintain the 
photoionization of hydrogen; it is lost from the system when 
the hydrogen recombines into levels above the ground state. 
By contrast, most of the energy of the stellar wind is con- 
verted into thermal energy of hot gas that accelerates the 
shell. The kinetic energy can be accumulated in the system 
for a long time unless it is dissipated, and the thermal energy 
of hot gas can also be saved if the density in the bubble is 
sufficiently low. For a plasma of solar chemical composition 
with T = 10 7 K and n — 10 -2 cm~ 3 , the cooling time in col- 
lisional ionization equilibrium is several times 10 7 yr. 
Although the LBV phase of the star induces very rapid 
changes in the energy balance of the circumstellar medium, 
its impact over longer timescales is limited as a result of the 
brevity of the LBV eruption. 

The above conclusions regarding the inefficiency of con- 
verting the energy flux of photoionizing UV photons into 
kinetic energy are modified in the presence of large-scale 
density gradients (i.e., at the edge of a molecular cloud). 
Generally speaking, the resulting “champagne flow” can 
lead to an efficiency of about 1% for converting the stellar 
UV flux into kinetic energy of expansion (see Yorke 1986). 
The exact value for efficiency depends on the details of the 
problem, however; it is conceivable in some cases that the 


champagne phase is negligibly short compared to the life- 
time of the star. We note that simulations that assess the role 
of stellar winds in the presence of champagne flows are cur- 
rently rare. Comeron (1997) examined the evolution of 
wind-driven H n regions across a density discontinuity for 
various parameters but without the consideration of energy 
transfer efficiencies. 

For the simulations of H n region evolution without a 
stellar wind the analytical prediction for the kinetic energy 
in the system at the end of the star’s lifetime differs from our 
numerical result by less than 30%. The analytical solution 
for the energy transfer efficiency in the case of the combined 
H ii region/SWB is just an upper limit and overestimates 
the transfer efficiencies into kinetic energy and thermal 
energy by factors of 6.7 and 9.7, respectively. 

This detailed examination is suitable to improve studies 
of the energization of the ISM in the solar neighborhood 
(e.g., Abbott 1982), which use energy transfer efficiencies as 
theoretical input. It is also important when considering the 
heating of the Galactic disk as a global phenomenon self- 
regulated by star formation. To address this, the effects of 
overlapping H n regions and SWBs in OB clusters and 
associations still need to be assessed, which is the subject of 
a future investigation. 

We thank Sabine Richling, Don Cox, Mordecai-Mark 
Mac Low, Jose Franco, Ralf-Jiirgen Dettmar, and Matthias 
Wrigge for helpful comments and interesting discussions. 
This work was supported by the Deutsche Forschungsge- 
meinschaft (DFG) under grant He 1487/17, by the 
Graduiertenkolleg GRK 118 in Bochum, and by the 
National Aeronautics and Space Administration (NASA) 
under grant NRA-99-01-ATP-065. Part of the research 
described in this paper was conducted at the Jet Propulsion 
Laboratory (JPL), California Institute of Technology. T. F. 
also gratefully acknowledges financial support by the 
German Academic Exchange Service (Doktorandenstipen- 
dium HSP II/AUFE) and the hospitality of the Department 
of Physics at the University of Wisconsin at Madison during 
a research visit. The computations were performed at the 
Rechenzentrum der Universitat Kiel, the Konrad- 
Zuse-Zentrum fur Informationstechnik in Berlin, and the 
John von Neumann-Institut fur Computing in Jiilich. We 
also wish to thank an anonymous referee for valuable 
comments. 


REFERENCES 


Abbott, D. C. 1982, ApJ, 263, 723 

Allen, C. W. 1973, Astrophysical Quantities (3d ed.; London; Athlone) 

Allen, D. A., & Burton, M. G. 1993, Nature, 363, 54 

Arthur, S. J„ Henney, W. J., & Dyson, J. E. 1996, A&A, 313, 897 

Avedisova, V. S. 1972, Soviet Astron., 15, 708 

Beltrametti, M., Tenorio-Tagle, G., & Yorke, H. W. 1982, A&A, 1 12, 1 

Benaglia, P., &Cappa, C. E. i 999, A&A, 346,979 

Bochkarev, N. G. 1988, Nature, 332, 518 

Brighenti, F., & D’Ercole, A. 1997, MNRAS, 285, 387 

Canto, J., Raga, A., Steffen, W., & Shapiro, P. R. 1998, ApJ, 502, 695 

Cappa, C. E., & Benaglia, P. 1998, AJ, 1 16, 1906 

Cappa, C. E., Niemela, V. S., Herbstmeier, U., & Koribalski, B. 1996, 
A&A, 3 12, 283 

Capriotti, E. R. 1973, ApJ, 179,495 

Castor, J., McCray, R., & Weaver, R. 1975, ApJ, 200, L107 

Chu, Y.-H. 1982, ApJ, 254, 578 

. 2000, Rev. Mexicana Astron. Astrofis. Ser. Conf., 9, 262 

Comeron, F. 1997, A&A, 326, 1 195 

Downes, D., Genzel, R., Becklin, E. E., & Wynn-Williams, C. G. 1981, 
ApJ, 244, 869 

Dyson, J. E. 1973, A&A, 23, 381 

Dyson, J. E., & de Vries, J. 1 972, A&A, 20, 223 


Falle, S. A. E. G. 1975, A&A, 43, 323 
Frank, A., Ryu, D., & Davidson, K. 1998, ApJ, 500, 291 
Garcia-Segura, G., & Franco, J. 1 996, ApJ, 469, 1 7 1 
Garcia-Segura, G., Langer, N., & Mac Low, M.-M. 1996a, A&A, 316, 133 
Garcia-Segura, G., & Mac Low, M.-M. 1995, ApJ, 455, 145 
Garcia-Segura, G., Mac Low, M.-M., & Langer, N. 1996b, A&A, 305, 229 
(GML1) 

Genzel, R., & Stutzki, J. 1989, ARA&A, 27, 41 
Hester, J. J.,etal. 1996,AJ, 111,2349 

Humphreys, R. M. 1991, in IAU Symp. 143, Wolf-Rayet Stars and Inter- 
relations with Other Massive Stars in Galaxies, ed. K. A. van der Hucht 
& B. Hidayat (Dordrecht: Kluwer), 485 
Koo, B.-C., & McKee, C. F. 1992a, ApJ, 388, 93 

. 1992b, ApJ, 388, 103 

Lasker, B. M. 1967, ApJ, 149, 23 
Levenson, N. A., et al. 2000, ApJ, 533, L53 
Lozinskaya, T. A. 1982, Ap&SS, 87, 313 

Mac Low, M.-M. 2000, Rev. Mexicana Astron. Astrofis. Ser. Conf., 9, 273 
Mac Low, M.-M., & Norman, M. L. 1993, ApJ, 407, 207 
Marston, A. P. 1996, AJ, 1 12, 2828 
- 1997, ApJ, 475, 188 

Martins, F., Schaerer, D., & Hillier, D. J. 2002, A&A, 382, 999 


910 


FREYER, HENSLER, & YORKE 


McCaughreari, M. J., & Mac Low, M.-M. 1997, AJ, 1 13, 391 
McKee, C. F., Van Buren, D., & Lazareff, B. 1984, ApJ, 278, LI 15 
Moffat, A. F. J., Firmani, C., McLean, L S., & Seggewiss, W. 1982, in 
IAU Symp. 99, Wolf-Rayet Stars: Observations, Physics, Evolution, ed. 
C. W. H. de Loore & A. J. Willis (Dordrecht: Reidel), 577 
Oey, M. S., & Massey, P. 1994, ApJ, 425, 635 

Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active 
Galactic Nuclei (Mill Valley: University Science Books) 

Pikel’ner, S. B. 1968, Astrophys. Lett., 2, 97 

Pilbratt, G. L., Altieri, B., Blommaert, J. A. D. L., Fridlund, C. V. M., 
Tauber, J. A., & Kessler, M. F. 1998, A&A, 333, L9 
Pittard, J. M., Dyson, J. E., & Hartquist, T. W. 2001a, A&A, 367, 1000 
Pittard, J. M., Hartquist, T. W., & Dyson, J. E. 2001b, A&A, 373, 1043 
Raymond, J. C., & Smith, B. W. 1977, ApJS, 35, 419 
Richling, S., & Yorke, H. W. 2000, ApJ, 539, 258 
Rozyczka, M. 1985, A&A, 143, 59 
Rozyczka, M., & Tenorio-Tagle, G. 1985a, A&A, 147, 202 

.1985b, A&A, 147,209 

Ryu, D., & Vishniac, E. T. 1988, ApJ, 331, 350 

Salas, L., et al. 1999, ApJ, 51 1, 822 

Sarazin, C. L., & White, R. E. 1987, ApJ, 320, 32 

Scoville, N., Kleinmann, S. G., Hall, D. N. B., & Ridgway, S. T. 1983, ApJ, 
275,201 


Soker, N. 1990,AJ,99, 1869 
— . 1998, MNRAS, 299, 562 

Spitzer, L. 1978, Physical Processes in the Interstellar Medium (New York: 
Wiley) 

Stone, J. M., Xu, J., & Mundy, L. G. 1995, Nature, 377, 315 
Strickland, D. K., & Stevens, I. R. 1998, MNRAS, 297, 747 
Taylor, K. N. R., Storey, J. W. V., Sandell, G., Williams, P. M., & Zealey, 
W. J. 1984, Nature, 31 1,236 

Tenorio-Tagle, G., Beltrametti, M., Bodenheimer, P., & Yorke, H. W. 
1982, A&A, 112, 104 

van Leer, B. 1977, J. Comput. Phys., 23, 276 
Vishniac, E. T. 1983, ApJ, 274, 152 
Vishniac, E. T., & Ryu, D. 1989, ApJ, 337, 917 

Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 
218, 377 (erratum 220, 742) 

Williams, R. J. R. 1999, MNRAS, 310, 789 
Wrigge, M. 1999, A&A, 343, 599 

Wrigge, M., Wendker, H. J., & Wisotzki, L. 1994, A&A, 286, 219 
Yang, H., Chu, Y.-H., Skillman, E. D., & Terlevich, R. 1996, AJ, 1 12, 146 
Yorke, H. W. 1 986, ARA&A, 24, 49 

Yorke, H. W., Bodenheimer, P., & Tenorio-Tagle, G. 1982, A&A, 108, 25 
Yorke, H. W., & Kaisig, M. 1995, Comput. Phys. Commun,, 89, 29 
Yorke, H. W., & Welz, A. 1996, A&A, 315, 555 



